home *** CD-ROM | disk | FTP | other *** search
/ Amiga Plus 1995 #5 & #6 / Amiga Plus CD - 1995 - No. 5 and 6.iso / pd / daten / astrolog / src / matrix.c < prev    next >
C/C++ Source or Header  |  1995-08-11  |  22KB  |  808 lines
  1. /*                                                               -*- C -*-
  2. ** Astrolog (Version 4.40) File: matrix.c
  3. **
  4. ** IMPORTANT NOTICE: The graphics database and chart display routines
  5. ** used in this program are Copyright (C) 1991-1995 by Walter D. Pullen
  6. ** (astara@u.washington.edu). Permission is granted to freely use and
  7. ** distribute these routines provided one doesn't sell, restrict, or
  8. ** profit from them in any way. Modification is allowed provided these
  9. ** notices remain with any altered or edited versions of the program.
  10. **
  11. ** The main planetary calculation routines used in this program have
  12. ** been Copyrighted and the core of this program is basically a
  13. ** conversion to C of the routines created by James Neely as listed in
  14. ** Michael Erlewine's 'Manual of Computer Programming for Astrologers',
  15. ** available from Matrix Software. The copyright gives us permission to
  16. ** use the routines for personal use but not to sell them or profit from
  17. ** them in any way.
  18. **
  19. ** The PostScript code within the core graphics routines are programmed
  20. ** and Copyright (C) 1992-1993 by Brian D. Willoughby
  21. ** (brianw@sounds.wa.com). Conditions are identical to those above.
  22. **
  23. ** The extended accurate ephemeris databases and formulas are from the
  24. ** calculation routines in the program "Placalc" and are programmed and
  25. ** Copyright (C) 1989,1991,1993 by Astrodienst AG and Alois Treindl
  26. ** (alois@azur.ch). The use of that source code is subject to
  27. ** regulations made by Astrodienst Zurich, and the code is not in the
  28. ** public domain. This copyright notice must not be changed or removed
  29. ** by any user of this program.
  30. **
  31. ** Initial programming 8/28,30, 9/10,13,16,20,23, 10/3,6,7, 11/7,10,21/1991.
  32. ** X Window graphics initially programmed 10/23-29/1991.
  33. ** PostScript graphics initially programmed 11/29-30/1992.
  34. ** Last code change made 1/29/1995.
  35. */
  36.  
  37. /* $VER: $Id: matrix.c,v 1.2 1995/07/02 22:21:59 tf Exp $ */
  38.  
  39. #include "astrolog.h"
  40.  
  41.  
  42. #ifdef MATRIX
  43. real MC, Asc, RA, OB;
  44.  
  45.  
  46. /*
  47. ******************************************************************************
  48. ** Assorted Calculations.
  49. ******************************************************************************
  50. */
  51.  
  52. /* Given a month, day, and year, convert it into a single Julian day value, */
  53. /* i.e. the number of days passed since a fixed reference date.             */
  54.  
  55. long MdyToJulian(mon, day, yea)
  56. int mon, day, yea;
  57. {
  58.  
  59. #ifndef PLACALC
  60.   long im, j;
  61.  
  62.   im = 12 * ( (long)yea+4800 ) + (long)mon-3;
  63.   j = ( 2*(im%12) + 7 + 365*im ) / 12;
  64.   j += (long)day + im/48 - 32083;
  65.  
  66.   if (j > 2299171)                   /* Take care of dates in */
  67.     j += im/4800 - im/1200 + 38;     /* Gregorian calendar.   */
  68.  
  69.   return j;
  70.  
  71. #else /* PLACALC */
  72.  
  73.   int fGreg = fTrue;
  74.  
  75.   if (yea < yeaJ2G || (yea == yeaJ2G && (mon < monJ2G || (mon == monJ2G && day < 15))))
  76.     fGreg = fFalse;
  77.  
  78.   return (long)RFloor(julday(mon, day, yea, 12.0, fGreg)+rRound);
  79.  
  80. #endif
  81. }
  82.  
  83.  
  84. /* Like above but return a fractional Julian time given the extra info. */
  85.  
  86. real MdytszToJulian(mon, day, yea, tim, dst, zon)
  87. int mon, day, yea;
  88. real tim, dst, zon;
  89. {
  90.   return (real)MdyToJulian(mon, day, yea) + (DecToDeg(tim) + DecToDeg(zon) - DecToDeg(dst)) / 24.0;
  91. }
  92.  
  93.  
  94. /* Take a Julian day value, and convert it back into the corresponding */
  95. /* month, day, and year.                                               */
  96.  
  97. void JulianToMdy(JD, mon, day, yea)
  98. real JD;
  99. int *mon, *day, *yea;
  100. {
  101.  
  102. #ifndef PLACALC
  103.   long L, N, IT, JT, K, IK;
  104.  
  105.   L  = (long)RFloor(JD+rRound)+68569L;
  106.   N  = Dvd(4L*L, 146097L);
  107.   L  -= Dvd(146097L*N + 3L, 4L);
  108.   IT = Dvd(4000L*(L+1L), 1461001L);
  109.   L  -= Dvd(1461L*IT, 4L) - 31L;
  110.   JT = Dvd(80L*L, 2447L);
  111.   K  = L-Dvd(2447L*JT, 80L);
  112.   L  = Dvd(JT, 11L);
  113.   JT += 2L - 12L*L;
  114.   IK = 100L*(N-49L) + IT + L;
  115.   *mon = (int)JT; *day = (int)K; *yea = (int)IK;
  116.  
  117. #else /* PLACALC */
  118.  
  119.   double tim;
  120.   revjul(JD, JD >= 2299171.0 /* October 15, 1582 */, mon, day, yea, &tim);
  121.  
  122. #endif
  123. }
  124.  
  125.  
  126. /* This is a subprocedure of CastChart(). Once we have the chart parameters, */
  127. /* calculate a few important things related to the date, i.e. the Greenwich  */
  128. /* time, the Julian day and fractional part of the day, the offset to the    */
  129. /* sidereal, and a couple of other things.                                   */
  130.  
  131. real ProcessInput(fDate)
  132. bool fDate;
  133. {
  134.   real Off, Ln;
  135.  
  136.   TT = RSgn(TT)*RFloor(RAbs(TT))+RFract(RAbs(TT))*100.0/60.0 + (DecToDeg(ZZ) - DecToDeg(SS));
  137.   OO = DecToDeg(OO);
  138.   AA = Min(AA, 89.9999);        /* Make sure the chart isn't being cast */
  139.   AA = Max(AA, -89.9999);       /* on the precise north or south pole.  */
  140.   AA = RFromD(DecToDeg(AA));
  141.  
  142.   /* if parameter 'fDate' isn't set, then we can assume that the true time */
  143.   /* has already been determined (as in a -rm switch time midpoint chart). */
  144.  
  145.   if (fDate) {
  146.     is.JD = (real)MdyToJulian(MM, DD, YY);
  147.  
  148.     if (!us.fProgress || us.fSolarArc)
  149.       T = (is.JD + TT/24.0 - 2415020.5) / 36525.0;
  150.  
  151.     else {
  152.       /* Determine actual time that a progressed chart is to be cast for. */
  153.  
  154.       T = ((is.JD + TT/24.0 + (is.JDp - (is.JD + TT/24.0)) / us.rProgDay) -
  155.         2415020.5) / 36525.0;
  156.     }
  157.   }
  158.  
  159.   /* Compute angle that the ecliptic is inclined to the Celestial Equator */
  160.   OB = RFromD(23.452294-0.0130125*T);
  161.  
  162.   Ln  = Mod((933060-6962911*T+7.5*T*T)/3600.0);    /* Mean lunar node */
  163.   Off = (259205536.0*T+2013816.0)/3600.0;          /* Mean Sun        */
  164.   Off = 17.23*RSin(RFromD(Ln))+1.27*RSin(RFromD(Off))-(5025.64+1.11*T)*T;
  165.   Off = (Off-84038.27)/3600.0;
  166.   is.rSid = (us.fSiderial ? Off : 0.0) + us.rZodiacOffset;
  167.   return Off;
  168. }
  169.  
  170.  
  171. /* Convert polar to rectangular coordinates. */
  172.  
  173. void PolToRec(A, R, X, Y)
  174. real A, R, *X, *Y;
  175. {
  176.   if (A == 0.0)
  177.     A = rSmall;
  178.  
  179.   *X = R*RCos(A);
  180.   *Y = R*RSin(A);
  181. }
  182.  
  183.  
  184. /* Convert rectangular to polar coordinates. */
  185.  
  186. void RecToPol(X, Y, A, R)
  187. real X, Y, *A, *R;
  188. {
  189.   if (Y == 0.0)
  190.     Y = rSmall;
  191.  
  192.   *R = RSqr(X*X + Y*Y);
  193.   *A = Angle(X, Y);
  194. }
  195.  
  196.  
  197. /* Convert rectangular to spherical coordinates. */
  198.  
  199. real RecToSph(B, L, O)
  200. real B, L, O;
  201. {
  202.   real R, Q, G, X, Y, A;
  203.  
  204.   A = B; R = 1.0;
  205.   PolToRec(A, R, &X, &Y);
  206.  
  207.   Q = Y; R = X; A = L;
  208.   PolToRec(A, R, &X, &Y);
  209.  
  210.   G = X; X = Y; Y = Q;
  211.   RecToPol(X, Y, &A, &R);
  212.  
  213.   A += O;
  214.   PolToRec(A, R, &X, &Y);
  215.  
  216.   Q = RAsin(Y);
  217.   Y = X; X = G;
  218.   RecToPol(X, Y, &A, &R);
  219.  
  220.   if (A < 0.0)
  221.     A += 2*rPi;
  222.  
  223.   G = A;
  224.   return G;  /* We only ever care about and return one of the coordinates. */
  225. }
  226.  
  227.  
  228. /* Do a coordinate transformation: Given a longitude and latitude value,    */
  229. /* return the new longitude and latitude values that the same location      */
  230. /* would have, were the equator tilted by a specified number of degrees.    */
  231. /* In other words, do a pole shift! This is used to convert among ecliptic, */
  232. /* equatorial, and local coordinates, each of which have zero declination   */
  233. /* in different planes. In other words, take into account the Earth's axis. */
  234.  
  235. void CoorXform(azi, alt, tilt)
  236. real *azi, *alt, tilt;
  237. {
  238.   real x, y, a1, l1;
  239.   real sinalt, cosalt, sinazi, sintilt, costilt;
  240.  
  241.   sinalt = RSin(*alt); cosalt = RCos(*alt); sinazi = RSin(*azi);
  242.   sintilt = RSin(tilt); costilt = RCos(tilt);
  243.  
  244.   x = cosalt * sinazi * costilt;
  245.   y = sinalt * sintilt;
  246.   x -= y;
  247.   a1 = cosalt;
  248.   y = cosalt * RCos(*azi);
  249.   l1 = Angle(y, x);
  250.   a1 = a1 * sinazi * sintilt + sinalt * costilt;
  251.   a1 = RAsin(a1);
  252.   *azi = l1; *alt = a1;
  253. }
  254.  
  255.  
  256. /* This is another subprocedure of CastChart(). Calculate a few variables */
  257. /* corresponding to the chart parameters that are used later on. The      */
  258. /* astrological vertex (object number twenty) is also calculated here.    */
  259.  
  260. void ComputeVariables(vtx)
  261. real *vtx;
  262. {
  263.   real R, R2, B, L, O, G, X, Y, A;
  264.  
  265.   RA = RFromD(Mod((6.6460656+2400.0513*T+2.58E-5*T*T+TT)*15.0-OO));
  266.   R2 = RA; O = -OB; B = AA; A = R2; R = 1.0;
  267.   PolToRec(A, R, &X, &Y);
  268.  
  269.   X *= RCos(O);
  270.   RecToPol(X, Y, &A, &R);
  271.  
  272.   MC = Mod(is.rSid + DFromR(A));              /* Midheaven */
  273.  
  274. #if FALSE
  275.   L = R2;
  276.   G = RecToSph(B, L, O);
  277.   Asc = Mod(is.rSid + Mod(G+rPiHalf));        /* Ascendant */
  278. #endif
  279.  
  280.   L= R2+rPi; B = rPiHalf-RAbs(B);
  281.  
  282.   if (AA < 0.0)
  283.     B = -B;
  284.  
  285.   G = RecToSph(B, L, O);
  286.   *vtx = Mod(is.rSid + DFromR(G+rPiHalf));    /* Vertex */
  287. }
  288.  
  289.  
  290. /*
  291. ******************************************************************************
  292. ** House Cusp Calculations.
  293. ******************************************************************************
  294. */
  295.  
  296. /* The following three functions calculate the Midheaven, Ascendant, and  */
  297. /* East Point of the chart in question, based on time and location. The   */
  298. /* first two are also used in some of the house cusp calculation routines */
  299. /* as a quick way to get the 10th and 1st cusps. The East Point object is */
  300. /* technically defined as the Ascendant's position at zero latitude.      */
  301.  
  302. real CuspMidheaven()
  303. {
  304.   real MC;
  305.  
  306.   MC = RAtn(RTan(RA)/RCos(OB));
  307.  
  308.   if (MC < 0.0)
  309.     MC += rPi;
  310.  
  311.   if (RA > rPi)
  312.     MC += rPi;
  313.  
  314.   return Mod(DFromR(MC)+is.rSid);
  315. }
  316.  
  317. real CuspAscendant()
  318. {
  319.   real Asc;
  320.  
  321.   Asc = Angle(-RSin(RA)*RCos(OB)-RTan(AA)*RSin(OB), RCos(RA));
  322.   return Mod(DFromR(Asc)+is.rSid);
  323. }
  324.  
  325. real CuspEastPoint()
  326. {
  327.   real EP;
  328.  
  329.   EP = Angle(-RSin(RA)*RCos(OB), RCos(RA));
  330.   return Mod(DFromR(EP)+is.rSid);
  331. }
  332.  
  333.  
  334. /* These are various different algorithms for calculating the house cusps: */
  335.  
  336. real CuspPlacidus(deg, FF, fNeg)
  337. real deg, FF;
  338. bool fNeg;
  339. {
  340.   real LO, R1, XS, X;
  341.   int i;
  342.  
  343.   R1 = RA+RFromD(deg);
  344.   X = fNeg ? 1.0 : -1.0;
  345.  
  346.   /* Looping 10 times is arbitrary, but it's what other programs do. */
  347.  
  348.   for (i = 1; i <= 10; i++) {
  349.  
  350.     /* This formula works except at 0 latitude (AA == 0.0). */
  351.  
  352.     XS = X*RSin(R1)*RTan(OB)*RTan(AA == 0.0 ? 0.0001 : AA);
  353.     XS = RAcos(XS);
  354.     if (XS < 0.0)
  355.       XS += rPi;
  356.     R1 = RA + (fNeg ? rPi-(XS/FF) : (XS/FF));
  357.   }
  358.  
  359.   LO = RAtn(RTan(R1)/RCos(OB));
  360.  
  361.   if (LO < 0.0)
  362.     LO += rPi;
  363.  
  364.   if (RSin(R1) < 0.0)
  365.     LO += rPi;
  366.  
  367.   return DFromR(LO);
  368. }
  369.  
  370.  
  371. void HousePlacidus()
  372. {
  373.   int i;
  374.  
  375.   house[1] = Mod(Asc-is.rSid);
  376.   house[4] = Mod(MC+rDegHalf-is.rSid);
  377.   house[5] = CuspPlacidus(30.0, 3.0, fFalse) + rDegHalf;
  378.   house[6] = CuspPlacidus(60.0, 1.5, fFalse) + rDegHalf;
  379.   house[2] = CuspPlacidus(120.0, 1.5, fTrue);
  380.   house[3] = CuspPlacidus(150.0, 3.0, fTrue);
  381.  
  382.   for (i = 1; i <= cSign; i++) 
  383.   {
  384.     if (i <= 6)
  385.       house[i] = Mod(house[i]+is.rSid);
  386.     else
  387.       house[i] = Mod(house[i-6]+rDegHalf);
  388.   }
  389. }
  390.  
  391. void HouseKoch()
  392. {
  393.   real A1, A2, A3, KN, D, X;
  394.   int i;
  395.  
  396.   A1 = RSin(RA)*RTan(AA)*RTan(OB);
  397.   A1 = RAsin(A1);
  398.  
  399.   for (i = 1; i <= cSign; i++) {
  400.     D = Mod(60.0+30.0*(real)i);
  401.     A2 = D/rDegQuad-1.0; KN = 1.0;
  402.     if (D >= rDegHalf) {
  403.       KN = -1.0;
  404.       A2 = D/rDegQuad-3.0;
  405.     }
  406.     A3 = RFromD(Mod(DFromR(RA)+D+A2*DFromR(A1)));
  407.     X = Angle(RCos(A3)*RCos(OB)-KN*RTan(AA)*RSin(OB), RSin(A3));
  408.     house[i] = Mod(DFromR(X)+is.rSid);
  409.   }
  410. }
  411.  
  412. void HouseEqual()
  413. {
  414.   int i;
  415.  
  416.   for (i = 1; i <= cSign; i++)
  417.     house[i] = Mod(Asc-30.0+30.0*(real)i);
  418. }
  419.  
  420. void HouseCampanus()
  421. {
  422.   real KO, DN, X;
  423.   int i;
  424.  
  425.   for (i = 1; i <= cSign; i++) 
  426.   {
  427.     KO = RFromD(60.000001+30.0*(real)i);
  428.     DN = RAtn(RTan(KO)*RCos(AA));
  429.  
  430.     if (DN < 0.0)
  431.       DN += rPi;
  432.  
  433.     if (RSin(KO) < 0.0)
  434.       DN += rPi;
  435.  
  436.     X = Angle(RCos(RA+DN)*RCos(OB)-RSin(DN)*RTan(AA)*RSin(OB), RSin(RA+DN));
  437.     house[i] = Mod(DFromR(X)+is.rSid);
  438.   }
  439. }
  440.  
  441. void HouseMeridian()
  442. {
  443.   real D, X;
  444.   int i;
  445.  
  446.   for (i = 1; i <= cSign; i++) 
  447.   {
  448.     D = RFromD(60.0+30.0*(real)i);
  449.     X = Angle(RCos(RA+D)*RCos(OB), RSin(RA+D));
  450.     house[i] = Mod(DFromR(X)+is.rSid);
  451.   }
  452. }
  453.  
  454. void HouseRegiomontanus()
  455. {
  456.   real D, X;
  457.   int i;
  458.  
  459.   for (i = 1; i <= cSign; i++) 
  460.   {
  461.     D = RFromD(60.0+30.0*(real)i);
  462.     X = Angle(RCos(RA+D)*RCos(OB)-RSin(D)*RTan(AA)*RSin(OB), RSin(RA+D));
  463.     house[i] = Mod(DFromR(X)+is.rSid);
  464.   }
  465. }
  466.  
  467. void HousePorphyry()
  468. {
  469.   real X, Y;
  470.   int i;
  471.  
  472.   X = Asc-MC;
  473.  
  474.   if (X < 0.0)
  475.     X += rDegMax;
  476.  
  477.   Y = X/3.0;
  478.  
  479.   for (i = 1; i <= 2; i++)
  480.     house[i+4] = Mod(rDegHalf+MC+i*Y);
  481.  
  482.   X = Mod(rDegHalf+MC)-Asc;
  483.  
  484.   if (X < 0.0)
  485.     X += rDegMax;
  486.  
  487.   house[1]=Asc;
  488.  
  489.   Y = X/3.0;
  490.  
  491.   for (i = 1; i <= 3; i++)
  492.     house[i+1] = Mod(Asc+i*Y);
  493.  
  494.   for (i = 1; i <= 6; i++)
  495.     house[i+6] = Mod(house[i]+rDegHalf);
  496. }
  497.  
  498. void HouseMorinus()
  499. {
  500.   real D, X;
  501.   int i;
  502.  
  503.   for (i = 1; i <= cSign; i++) 
  504.   {
  505.     D = RFromD(60.0+30.0*(real)i);
  506.     X = Angle(RCos(RA+D), RSin(RA+D)*RCos(OB));
  507.     house[i] = Mod(DFromR(X)+is.rSid);
  508.   }
  509. }
  510.  
  511. real CuspTopocentric(deg)
  512. real deg;
  513. {
  514.   real OA, X, LO;
  515.  
  516.   OA = ModRad(RA+RFromD(deg));
  517.   X = RAtn(RTan(AA)/RCos(OA));
  518.   LO = RAtn(RCos(X)*RTan(OA)/RCos(X+OB));
  519.  
  520.   if (LO < 0.0)
  521.     LO += rPi;
  522.  
  523.   if (RSin(OA) < 0.0)
  524.     LO += rPi;
  525.  
  526.   return LO;
  527. }
  528.  
  529. void HouseTopocentric()
  530. {
  531.   real TL, P1, P2, LT;
  532.   int i;
  533.  
  534.   house[4] = ModRad(RFromD(MC+rDegHalf-is.rSid));
  535.   TL = RTan(AA); P1 = RAtn(TL/3.0); P2 = RAtn(TL/1.5); LT = AA;
  536.   AA = P1; house[5] = CuspTopocentric(30.0) + rPi;
  537.   AA = P2; house[6] = CuspTopocentric(60.0) + rPi;
  538.   AA = LT; house[1] = CuspTopocentric(90.0);
  539.   AA = P2; house[2] = CuspTopocentric(120.0);
  540.   AA = P1; house[3] = CuspTopocentric(150.0);
  541.   AA = LT;
  542.   for (i = 1; i <= 6; i++) 
  543.   {
  544.     house[i] = Mod(DFromR(house[i])+is.rSid);
  545.     house[i+6] = Mod(house[i]+rDegHalf);
  546.   }
  547. }
  548.  
  549.  
  550. /*
  551. ******************************************************************************
  552. ** Planetary Position Calculations.
  553. ******************************************************************************
  554. */
  555.  
  556. /* Given three values, return them combined as the coefficients of a */
  557. /* quadratic equation as a function of the chart time.               */
  558.  
  559. real ReadThree(r0, r1, r2)
  560. real r0, r1, r2;
  561. {
  562.   return RFromD(r0 + r1*T + r2*T*T);
  563. }
  564.  
  565.  
  566. /* Another coordinate transformation. This is used by the ComputePlanets() */
  567. /* procedure to rotate rectangular coordinates by a certain amount.        */
  568.  
  569. void RecToSph2(AP, AN, IN, X, Y, G)
  570. real AP, AN, IN, *X, *Y, *G;
  571. {
  572.   real R, D, A;
  573.  
  574.   RecToPol(*X, *Y, &A, &R); A += AP; PolToRec(A, R, X, Y);
  575.   D = *X; *X = *Y; *Y = 0.0; RecToPol(*X, *Y, &A, &R);
  576.   A += IN; PolToRec(A, R, X, Y);
  577.   *G = *Y; *Y = *X; *X = D; RecToPol(*X, *Y, &A, &R); A += AN;
  578.  
  579.   if (A < 0.0)
  580.     A += 2.0*rPi;
  581.  
  582.   PolToRec(A, R, X, Y);
  583. }
  584.  
  585.  
  586. /* Calculate some harmonic delta error correction factors to add onto the */
  587. /* coordinates of Jupiter through Pluto, for better accuracy.             */
  588.  
  589. void ErrorCorrect(ind, x, y, z)
  590. int ind;
  591. real *x, *y, *z;
  592. {
  593.   real U, V, W, A, S0, T0[4], FAR *pr;
  594.   int IK, IJ, irError;
  595.  
  596.   irError = errorcount[ind-oJup];
  597.  
  598.   pr = (lpreal)&errordata[erroroffset[ind-oJup]];
  599.  
  600.   for (IK = 1; IK <= 3; IK++) {
  601.     if (ind == oJup && IK == 3) {
  602.       T0[3] = 0.0;
  603.       break;
  604.     }
  605.     if (IK == 3)
  606.       irError--;
  607.     S0 = ReadThree(pr[0], pr[1], pr[2]); pr += 3;
  608.     A = 0.0;
  609.     for (IJ = 1; IJ <= irError; IJ++) {
  610.       U = *pr++; V = *pr++; W = *pr++;
  611.       A += RFromD(U)*RCos((V*T+W)*rPi/rDegHalf);
  612.     }
  613.     T0[IK] = DFromR(S0+A);
  614.   }
  615.   *x += T0[2]; *y += T0[1]; *z += T0[3];
  616. }
  617.  
  618.  
  619. /* Another subprocedure of the ComputePlanets() routine. Convert the final */
  620. /* rectangular coordinates of a planet to zodiac position and declination. */
  621.  
  622. void ProcessPlanet(ind, aber)
  623. int ind;
  624. real aber;
  625. {
  626.   real ang, rad;
  627.  
  628.   RecToPol(spacex[ind], spacey[ind], &ang, &rad);
  629.   planet[ind] = Mod(DFromR(ang) /*+ NU*/ - aber + is.rSid);
  630.   RecToPol(rad, spacez[ind], &ang, &rad);
  631.  
  632.   if (us.objCenter == oSun && ind == oSun)
  633.     ang = 0.0;
  634.  
  635.   ang = DFromR(ang);
  636.  
  637.   while (ang > rDegQuad)    /* Ensure declination is from -90..+90 degrees. */
  638.     ang -= rDegHalf;
  639.  
  640.   while (ang < -rDegQuad)
  641.     ang += rDegHalf;
  642.  
  643.   planetalt[ind] = ang;
  644. }
  645.  
  646.  
  647. /* This is probably the heart of the whole program of Astrolog. Calculate  */
  648. /* the position of each body that orbits the Sun. A heliocentric chart is  */
  649. /* most natural; extra calculation is needed to have other central bodies. */
  650.  
  651. void ComputePlanets()
  652. {
  653.   real helio[oNorm+1], 
  654.        helioalt[oNorm+1], 
  655.        helioret[oNorm+1],
  656.        heliox[oNorm+1], 
  657.        helioy[oNorm+1], 
  658.        helioz[oNorm+1];
  659.  
  660.   real aber = 0.0, AU, E, EA, E1, M, XW, YW, AP, AN, IN, X, Y, G, XS, YS, ZS;
  661.  
  662.   int ind = oSun, i;
  663.   OE *poe;
  664.  
  665.   while (ind <= (us.fUranian ? oNorm : cPlanet)) {
  666.  
  667.     if (ignore[ind] && ind > oSun)
  668.       goto LNextPlanet;
  669.  
  670.     poe = &rgoe[IoeFromObj(ind)];
  671.  
  672.     EA = M = ModRad(ReadThree(poe->ma0, poe->ma1, poe->ma2));
  673.     E = DFromR(ReadThree(poe->ec0, poe->ec1, poe->ec2));
  674.     for (i = 1; i <= 5; i++)
  675.       EA = M+E*RSin(EA);            /* Solve Kepler's equation */
  676.     AU = poe->sma;                  /* Semi-major axis         */
  677.     E1 = 0.01720209/(pow(AU, 1.5)*
  678.       (1.0-E*RCos(EA)));                     /* Begin velocity coordinates */
  679.     XW = -AU*E1*RSin(EA);                    /* Perifocal coordinates      */
  680.     YW = AU*E1*pow(1.0-E*E,0.5)*RCos(EA);
  681.     AP = ReadThree(poe->ap0, poe->ap1, poe->ap2);
  682.     AN = ReadThree(poe->an0, poe->an1, poe->an2);
  683.     IN = ReadThree(poe->in0, poe->in1, poe->in2); /* Calculate inclination  */
  684.     X = XW; Y = YW;
  685.     RecToSph2(AP, AN, IN, &X, &Y, &G);  /* Rotate velocity coords */
  686.     heliox[ind] = X; helioy[ind] = Y;
  687.     helioz[ind] = G;                    /* Helio ecliptic rectangtular */
  688.     X = AU*(RCos(EA)-E);                /* Perifocal coordinates for        */
  689.     Y = AU*RSin(EA)*pow(1.0-E*E,0.5);   /* rectangular position coordinates */
  690.     RecToSph2(AP, AN, IN, &X, &Y, &G);  /* Rotate for rectangular */
  691.     XS = X; YS = Y; ZS = G;             /* position coordinates   */
  692.     if (FBetween(ind, oJup, oPlu))
  693.       ErrorCorrect(ind, &XS, &YS, &ZS);
  694.     ret[ind] =                                        /* Helio daily motion */
  695.       (XS*helioy[ind]-YS*heliox[ind])/(XS*XS+YS*YS);
  696.     spacex[ind] = XS; spacey[ind] = YS; spacez[ind] = ZS;
  697.     ProcessPlanet(ind, 0.0);
  698. LNextPlanet:
  699.     ind += (ind == oSun ? 2 : (ind != cPlanet ? 1 : uranLo-cPlanet));
  700.   }
  701.   spacex[0] = spacey[0] = spacez[0] = 0.0;
  702.  
  703.   if (!us.objCenter) {
  704.     if (us.fVelocity)
  705.       for (i = 0; i <= oNorm; i++)    /* Use relative velocity */
  706.         ret[i] = RFromD(1.0);         /* if -v0 is in effect.  */
  707.     return;
  708.   }
  709. #endif /* MATRIX */
  710.  
  711.   /* A second loop is needed for geocentric charts or central bodies other */
  712.   /* than the Sun. For example, we can't find the position of Mercury in   */
  713.   /* relation to Pluto until we know the position of Pluto in relation to  */
  714.   /* the Sun, and since Mercury is calculated first, another pass needed.  */
  715.  
  716.   ind = us.objCenter;
  717.  
  718.   for (i = 0; i <= oNorm; i++) {
  719.     helio[i]    = planet[i];
  720.     helioalt[i] = planetalt[i];
  721.     helioret[i] = ret[i];
  722.     if (i != oMoo && i != ind) {
  723.       spacex[i] -= spacex[ind];
  724.       spacey[i] -= spacey[ind];
  725.       spacez[i] -= spacez[ind];
  726.     }
  727.   }
  728.   spacex[ind] = spacey[ind] = spacez[ind] = 0.0;
  729.   SwapR(&spacex[0], &spacex[ind]);
  730.   SwapR(&spacey[0], &spacey[ind]);    /* Do some swapping - we want   */
  731.   SwapR(&spacez[0], &spacez[ind]);    /* the central body to be in    */
  732.   SwapR(&spacex[1], &spacex[ind]);    /* object position number zero. */
  733.   SwapR(&spacey[1], &spacey[ind]);
  734.   SwapR(&spacez[1], &spacez[ind]);
  735.  
  736.   for (i = 1; i <= (us.fUranian ? oNorm : cPlanet);
  737.       i += (i == 1 ? 2 : (i != cPlanet ? 1 : uranLo-cPlanet))) {
  738.     if (ignore[i] && i > oSun)
  739.       continue;
  740.     XS = spacex[i]; YS = spacey[i]; ZS = spacez[i];
  741.     if (ind != oSun || i != oSun)
  742.       ret[i] = (XS*(helioy[i]-helioy[ind])-YS*(heliox[i]-heliox[ind]))/
  743.         (XS*XS+YS*YS);
  744.     if (ind == oSun)
  745.       aber = 0.0057756*RSqr(XS*XS+YS*YS+ZS*ZS)*DFromR(ret[i]); /* Aberration */
  746.     ProcessPlanet(i, aber);
  747.     if (us.fVelocity)                         /* Use relative velocity */
  748.       ret[i] = RFromD(ret[i]/helioret[i]);    /* if -v0 is in effect   */
  749.   }
  750. }
  751.  
  752.  
  753. #ifdef MATRIX
  754. /*
  755. ******************************************************************************
  756. ** Lunar Position Calculations
  757. ******************************************************************************
  758. */
  759.  
  760. /* Calculate the position and declination of the Moon, and the Moon's North  */
  761. /* Node. This has to be done separately from the other planets, because they */
  762. /* all orbit the Sun, while the Moon orbits the Earth.                       */
  763.  
  764. void ComputeLunar(moonlo, moonla, nodelo, nodela)
  765. real *moonlo, *moonla, *nodelo, *nodela;
  766. {
  767.   real LL, G, N, G1, D, L, ML, L1, MB, T1, Y, M = 3600.0;
  768.  
  769.   LL = 973563.0+1732564379.0*T-4.0*T*T;  /* Compute mean lunar longitude     */
  770.   G = 1012395.0+6189.0*T;                /* Sun's mean longitude of perigee  */
  771.   N = 933060.0-6962911.0*T+7.5*T*T;      /* Compute mean lunar node          */
  772.   G1 = 1203586.0+14648523.0*T-37.0*T*T;  /* Mean longitude of lunar perigee  */
  773.   D = 1262655.0+1602961611.0*T-5.0*T*T;  /* Mean elongation of Moon from Sun */
  774.   L = (LL-G1)/M; L1 = ((LL-D)-G)/M;      /* Some auxiliary angles            */
  775.   T1 = (LL-N)/M; D = D/M; Y = 2.0*D;
  776.  
  777.   /* Compute Moon's perturbations. */
  778.  
  779.   ML = 22639.6*RSinD(L) - 4586.4*RSinD(L-Y) + 2369.9*RSinD(Y) +
  780.     769.0*RSinD(2.0*L) - 669.0*RSinD(L1) - 411.6*RSinD(2.0*T1) -
  781.     212.0*RSinD(2.0*L-Y) - 206.0*RSinD(L+L1-Y);
  782.   ML += 192.0*RSinD(L+Y) - 165.0*RSinD(L1-Y) + 148.0*RSinD(L-L1) -
  783.     125.0*RSinD(D) - 110.0*RSinD(L+L1) - 55.0*RSinD(2.0*T1-Y) -
  784.     45.0*RSinD(L+2.0*T1) + 40.0*RSinD(L-2.0*T1);
  785.  
  786.   *moonlo = G = Mod((LL+ML)/M+is.rSid);    /* Lunar longitude */
  787.  
  788.   /* Compute lunar latitude. */
  789.  
  790.   MB = 18461.5*RSinD(T1) + 1010.0*RSinD(L+T1) - 999.0*RSinD(T1-L) -
  791.     624.0*RSinD(T1-Y) + 199.0*RSinD(T1+Y-L) - 167.0*RSinD(L+T1-Y);
  792.   MB += 117.0*RSinD(T1+Y) + 62.0*RSinD(2.0*L+T1) -
  793.     33.0*RSinD(T1-Y-L) - 32.0*RSinD(T1-2.0*L) - 30.0*RSinD(L1+T1-Y);
  794.   *moonla = MB =
  795.     RSgn(MB)*((RAbs(MB)/M)/rDegMax-RFloor((RAbs(MB)/M)/rDegMax))*rDegMax;
  796.  
  797.   /* Compute position of the North Lunar Node, either True or Mean. */
  798.  
  799.   if (us.fTrueNode)
  800.     N = N+5392.0*RSinD(2.0*T1-Y)-541.0*RSinD(L1)-442.0*RSinD(Y)+
  801.       423.0*RSinD(2.0*T1)-291.0*RSinD(2.0*L-2.0*T1);
  802.   *nodelo = Mod(N/M+is.rSid);
  803.   *nodela = 0.0;
  804. }
  805. #endif /* MATRIX */
  806.  
  807. /* matrix.c */
  808.